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, We propose explicit symplectic integrators of molecular dynamics (MD) algorithms for rigid- 

. body molecules in the canonical and isothermal-isobaric ensembles. We also present a symplectic 

algorithm in the constant normal pressure and lateral surface area ensemble and that combined with 
the Parrinello-Rahman algorithm. Employing the symplectic integrators for MD algorithms, there 
is a conserved quantity which is close to Hamiltonian. Therefore, we can perform a MD simulation 
• more stably than by conventional nonsymplectic algorithms. We applied this algorithm to a TIP3P 

pure water system at 300 K and compared the time evolution of the Hamiltonian with those by the 
nonsymplectic algorithms. We found that the Hamiltonian was conserved well by the symplectic 
algorithm even for a time step of 4 fs. This time step is longer than typical values of 0.5-2 fs which 
are used by the conventional nonsymplectic algorithms. 
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There are two models for molecules in molecular dynamics (MD) simulations. One model is a rigid-body model 
and the other is a flexible model. Relative coordinates in a molecule are fixed in the rigid-body model, while they 
vary in the flexible model. Because degrees of freedom in the rigid-body model are fewer than in the flexible model, 
the simulational cost is less expensive. Several MD techniques have thus been proposed for rigid-body molecules. 

One possibility for the rigid-body modeling is to constrain a bond length and a bond angle among atoms in the 
molecules such as in the SHAKE algorithm [lj . Although it is easy to write a computational program for this constraint 
algorithm, it requires iteration procedures to fulfill the constraint. It means that one has to perform implicit time 
development. 

' Another algorithm is a quaternion scheme which gives explicit time development. One integrator to carry out a 
quaternion MD simulation is the Gear's predictor-corrector algorithm 0. However, this algorithm is not a symplectic 
integrator Q nor time reversible. It hardly reflects characteristics of Hamiltonian dynamics. Another algorithm for 
the quaternion MD was proposed by Matubayasi and Nakahara Q|. Although this algorithm is not a symplectic 
integrator, it conserves volume in phase space and is time reversible. Miller et al. recently proposed a symplectic 
quaternion algorithm p|. This algorithm also conserves volume in phase space and is time reversible. However, this 
symplectic quaternion algorithm has been proposed only in the microcanonical ensemble. There is no symplectic 
quaternion algorithm in the canonical ensemble and in the isothermal-isobaric ensemble. 

A representative MD algorithm to obtain the canonical ensemble is the Nose thermostat 0, Q • Because the original 
Nose Hamiltonian gives dynamics in virtual time, a symplectic canonical MD simulation can be carried out in virtual 
time. However, a symplectic MD simulation cannot be realized in real time. Nonsymplectic integrators such as Gear's 
O ■ predictor-corrector algorithm are often employed for real-time development for the Nose thermostat. Hoover improved 
the Nose thermostat to propose the Nose- Hoover thermostat 0- Because the Nose-Hoover thermostat is not based 
on a Hamiltonian, there is no symplectic algorithm [t| for the Nose-Hoover thermostat. However, there exists an 
explicit time reversible integrator, although it does not conserve the volume in the phase space. This integrator was 
proposed by Martyna et al. [Toj |. Bond et al. then proposed a symplectic constant temperature algorithm in real 
time, which is refereed to as the Nose-Poincare thermostat ^lj. However, the original symplectic algorithm for the 
Nose-Poincare thermostat is an implicit integrator. Iterations are necessary for the thermostat. Nose improved the 
original algorithm and proposed an explicit symplectic integrator for the Nose-Poincare thermostat |l2 ] . Although 
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this formalism may not be widely known, we found it very powerful and useful as was shown in Refs. |13l Il4j and will 
be demonstrated below. 

An explicit MD algorithm closest to a svmplectic algorithm for rigid-body molecules in the canonical ensemble 
proposed so far is a combined algorithm |l5| of the svmplectic quaternion algorithm by Miller et al. |5j and time 
reversible algorithm for the Nose-Hoover thermostat Because the Nose-Hoover thermostat is a nonsymplectic 

algorithm, the whole algorithm is also nonsymplectic. Employing a symplectic MD algorithm, there is a conserved 
quantity which is close to Hamiltonian and the long-time deviation of the Hamiltonian is suppressed. Symplectic MD 
algorithms are thus getting popular recently. 

In this article, we propose an explicit symplectic MD algorithm for rigid-body molecules in the canonical ensemble 
by combining the quaternion algorithm by Miller et al. |j| with the explicit symplectic algorithm for the Nose-Poincare 
thermostat by Nose 01 • We further combine our algorithm with the Andersen barostat 0] to present an explicit 
symplectic MD algorithm for rigid-body molecules in the isothermal-isobaric ensemble. An explicit s ymp lectic MD 
algorithm for spherical atoms in the isothermal-isobaric ensemble has been presented in references [13, uM- The 
isothermal-isobaric algorithm in this article is an extension of the algorithm for spherical atoms to that for rigid-body 
molecules. We also present a symplectic integrator in the constant normal pressure and lateral surface area ensemble 
and a symplectic integrator combined with the Parrinello-Rahman algorithm. 

In Section [H] we first give brief reviews of the Nose-Poincare thermostat and the rigid-body MD algorithm. We 
then explain the symplectic MD algorithms for rigid-body molecules in the canonical, isothermal-isobaric, and related 
ensembles. In Section IIIII we compare our symplectic MD algorithm with nonsymplectic MD algorithms in the 
canonical ensemble. We apply our symplectic MD algorithm to a rigid-body water model and make numerical 
comparisons with the nonsymplectic MD algorithms. Section HVI is devoted to conclusions. 



II. METHODS 



Nose-Poincare thermostat 



The Nose-Poincare Hamiltonian -Hnp for N spherical atoms at temperature To is given by 0, 0] 

N -„/2 p2 

E(r^) + ^ + gk B T log s - H 
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where p\ and P s are the conjugate momenta for the coordinate of particle i and Nose's additional degree of 
freedom s, respectively. We have introduced a simplified notation by the superscript {N} for the set of coordinate 
and momentum vectors: r^ N > = (ri,r 2 , • • • ,t*jv) T and p?' N } = (p' l7 p2, ■ ■ ■ ,p' N ) T , where the superscript T stands 
for transpose. The real momentum p t and the virtual momentum p[ are related by 



Pi = Pi/s. 



(2) 



E is the potential energy. The constant rrii is the mass of particle i and Q is the artificial "mass" associated with s. 
The constant g corresponds to the number of degrees of freedom. In the case of a spherical atomic system, g equals 
3N (g equals 6N in the case of a rigid-body molecular system) . The Hamiltonian i?N is the original Nose Hamiltonian 
and Hq is the initial value of i?N- 

The equations of motion for the Nose-Poincare thermostat are given by 

Pi 
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where the dot above each variable stands for the time derivative and the relation of 

H N -H Q = (7) 
is used because H^ is conserved. Equations Q-© are the same as those for the Nose thermostat in the real time. 
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B. Molecular dynamics algorithm for rigid-body molecules in the microcanonical ensemble 

Hamiltonian for rigid-body molecules Hb.b are given by 0, 



N 



(8) 



where q i is a quaternion of molecule i, which indicates the orientation of the rigid-body molecule. Here, the quaternion 
q = (qo, <7i, q 2 , 93) T is related to the Euler angle (</), 6, tp) as follows: 



Qo = cos ( - ] cos 



qi = sin I - | cos 



q 2 = sin [ - | sin 



93 = cos \ - ] sin 
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The elements of the matrix S(q) are given by 

S(q) 
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The variable iti is the conjugate momentum for q i . The matrix D is a 4x4 matrix consisting of the inverse of the 
principal moments of inertia I\, I 2 , and I3 of molecule i: 



D 
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V / 3 "V 



(14) 



where Jo is an artificial constant. Note that the correct equations of motion for rigid-body molecules are obtained in 
the limit of Iq —>■ 00. In order to write the equations of motion more elegantly, we may introduce the angular velocity 



and the four-dimensional angular velocity 



U3 = (wi,W 2 ,W 3 ) T , 

> (4) = (0, uii, uj 2 , w 3 ) T 



(15) 
(16) 



where cj±, uj 2 , and o->3 are the angular velocities along each of the corresponding principal axes. In the limit of Iq — » 00, 
the four-dimensional angular velocity u> is related to 7Tj by 



.(4) 



1<- 
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In this limit the equations of motion for rigid-body molecules are obtained as follows: 



1 



S(qM 



(4) 



(17) 

(18) 
(19) 

Equation i|19|) is called the Euler equation of motion. Here, J is the 3x3 diagonal matrix whose diagonal elements 
are I\, I 2 , and I3. The vector 7V^ is the torque of molecule i, which is calculated by 



IiUJi = Ni - Wj X IiU 



Ni =^2r a x F a , 



(20) 
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where F a and r a are the coordinate and force of atom a, respectively, in a rigid-body- fixed coordinate system for 
molecule i. The torque TV^ is related to the potential energy E by 



where 




C. Symplectic molecular dynamics algorithm for rigid-body molecules combined with the Nose-Poincare 

thermostat 



We here present the explicit symplectic MD algorithm for rigid-body molecules in the canonical ensemble. We 
combine the Nose-Poincare Hamiltonian in Eq- CD El E2 

and the Hamiltonian for rigid-body molecules in Eq. (JSJ 
0, Il7| . The Nose-Poincare Hamiltonian for rigid-body molecules is given by 



#np-rb — s 



' N ,2 N ^ T 
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where r^ N ^ = (n, r"2, • • • , J"tv) t stands for the set of the coordinates of the center of mass for the rigid-body molecules. 
The vector tt[ is the conjugate momentum for quaternion q,. The real momentum 7r^ of the quaternion is related to 
the virtual momentum ir^ by 



s 



The equations of motion are given from the Hamiltonian in Eq. (|23() by 



Vi = 

Pi = 
Qi = 



rrii 



F t ~ -Pi , 

s 



(4) 



I l U) l = Ni- U>i X yliUij - -li^i , 
P, 



Q 
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The time development of a physical quantity Z(T) in the phase space T = (r^ N \p'^ N \ q^ N \n'^ N \ s, P S ) T is 
written by 



dZ _■ dZ 
~dt ~ ' df ' 

The formal solution of the time development of Z from time t to t + At is given by 

Z(t + At) = e DAt Z(t) , 

where e DAt is called a time propagator. The operator D is defined by 

8 



D = T ■ 



dT 



(31) 
(32) 

(33) 
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In the symplectic algorithm, the Hamiltonian in Eq. (12;SI) is separated into six terms here as follows: 



where 



#NP-RB — -^NP-RBO + -HnP-RBI + ^NP-RB2 
+ -f^NP — RB3 + -ffNP-RB4 + #NP-RB5 
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V%q = (-92,-93, 9o, 9i) T , 

^39 = (-93, 92,-9i, 9o) T - 



5/c B T log s - H 
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In the limit of Jo ~* °°, ^np-rbo goes to zero: .Hnp-rbo — * 0. Hereafter, then only Hamiltonians from -//np-rbi to 
i?NP-RB5 are considered. The second-order formula with respect to At is obtained by the decomposition of the time 
propagator exp [DAt] into a product of five time propagators: 
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Higher-order formulae can also be obtained in a similar manner. The explicit form of each operator is as follows: 
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There is no term higher than the second power of At in the time developments by -D4, because there is no conjugate 
pair in Pnp-rb4- Although there is a conjugate pair of and 7rJ in Pnp-rbi, the time developments of tj^ and tv[ 
by .Hnp-rbi are given by @ 



exp [DiAt] q, = cos (&iAt) q t + sin (&i At) Viqi , 
exp [£>i At] tx\ = cos (Cu At) tt- + sin (Oi At) Pi< , 



where 



47 lS 
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The time developments of q { and 71^ by P/2 and D3 are also obtained in the same way. Although there is another 
conjugate pair of s and P s in Pnp-rbs, the time developments of s and P s by D5 are given explicitly by [H 



• •x}) [P.-, A'] ••< = »•( J + 

<'X 1 )[./l;,A/]/b : : ( 1 + ^At 



(54) 
(55) 
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Finally, the explicit symplectic time developments for rigid-body molecules in the canonical ensemble is obtained 
from Eq. I|45|l . Here, a symbol of <— stands for a substitution in a computational program (i.e., the variables in each 
step adopt the substitutions in the preceding steps): 
Step 1. exp[D^At/2] operation: 



Step 2. explD^At/2] operation: 



Step 3. exp[D^At/2] operation: 



Step 4. exp[Z)2A£/2] operation: 



Step 5. exp[Z?iAi] operation: 
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Step 6. exp[£>2Ai/2] operation: 



Step 7. exp[D 3 Ai/2] operation: 



Step 8. exp[Z?4At/2] operation: 



Step 9. exp [D5A£/2] operation: 
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D. Symplectic molecular dynamics algorithm for rigid-body molecules combined with the Nose-Poincare 

thermostat and the Andersen barostat 



In this subsection we present the explicit symplectic MD algorithm for rigid-body molecules in the isothermal- 
isobaric ensemble. Hamiltonian for rigid-body molecules at temperature To and pressure Pq is given by combing the 
Hamiltonian in Eq. I|23|l and the Andersen barostat 0] as follows: 
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where p i and f i are the scaled momentum and the scaled coordinate by volume V and the degree of the Nose-Poincare 
thermostat s. They are related to p t and r, ; by 
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The constant W is the "mass" associated with V. The variable Py is the conjugate momenta for V. The constant 
H here is the initial value of the Nose- Andersen Hamiltonian H^a- The equations of motion are given by 
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where the relation of 



-Hna — Hq — 

is used. 

The Hamiltonian in the isothermal-isobaric ensemble is separated into six terms as follows: 

^NPA-RB = -HnPA-RBI + -HnPA-RB2 + #NPA-RB3 
+ #NPA-RB4 + -#NPA-RB5 + -HnPA— RB6 i 
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where the term of sJ^iLi [i^^'PaQi ) /8Ias 2 has been neglected again because it is zero in the limit of Iq — ► oo. 
As in the decomposition in Eq. 1|45|) in the canonical ensemble, the second-order formula is obtained for the time 
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propagator exp [D At] as a product of six time propagators: 

exp [DAt] = exp D 6 
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where Di, Z?2, ■ ■ ■ ^ are the time propagators which correspond to i?NPA— rbi j -f^NPA— RB2: ■ • ■ , Pnpa— RB61 respec- 
tively. 

According to the decomposition in Eq. i|10(j|) . the explicit symplectic time developments for rigid-body molecules 
in the isothermal-isobaric ensemble are given as follows: 
Step 1. exp[_DgAt/2] operation: 



Step 2. exp[D$At/2] operation: 



Step 3. exp[£>4At/2] operation: 



Step 4. exp[Z? 3 At/2] operation: 
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Step 5. exp[D 2 At/2] operation: 
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<l l <- cos f Ci2^ ) q L + sin f ( t 2^ J P 2 g, , (120) 



2 y 2 

^ - cos ( Ci 2^J ttJ + sin Ui2^-) $2< , (121) 



Step 6. exp[DiAt] operation: 
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/ N „ 2 N \ 

P s ^ P s + V - P, „ a + V - <7fc B T log a + H - gk B T a At , (127) 

AT _ 2 
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Step 7. exp[_D 2 Ai/2] operation: 



C,2 - ^-Trf V 2 q t , (129) 
<?* <- cos Ki2^ ) <?; + sin K»2^ ) , (130) 
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Step 8. exp[_D 3 Ai/2] operation: 



Step 9. exp[_D4At/2] operation: 



P s - P.+ ^^jl. (132) 
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Step 10. exp[Z?5At/2] operation: 



Step 11. exp[D 6 At/2] operation: 
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E. Symplectic molecular dynamics algorithm for rigid-body molecules in the constant temperature, 
constant normal pressure, and constant lateral surface area ensemble 



An explicit symplectic MD algorithm in the constant temperature, constant normal pressure, and constant lateral 
surface area ensemble is also easily obtained. In section lllDl the Andersen's constant pressure algorithm was employed 
for all three side lengths of the simulation cell. On the other hand, one of the side lengths of the simulation cell 
fluctuates in the constant normal pressure and constant lateral surface area ensemble. This ensemble is often used 
for membrane systems |l5(. The Hamiltonian for this ensemble is given by 



Pnpai-rb 



Pxi 



' N 

y „ , / 

' 2rriiS 2 L 2 ^— ' 2rriiS 2 

,i—l i—1 



N 



i=l 



+ BfiW,^,^,^,!) + g + gk B T logs + || + P Q AL H 



(145) 



where the variable Pl is the conjugate momenta for the side length L of the simulation cell along the x-axis. The 
constant A is the lateral surface area on the yz-plane. Therefore the volume of the simulation cell V is given by 
AL. Note that x components of l p i and are scaled by p X i — p X i/sL and Xi — Lxi , respectively, whereas y and z 
components of p i are scaled by Eq. J5J). The equations of motion are given by 



Pxi 
Pyi 



Pxi 

rrii 
Pyi 
m l 

F 



L 



Pzi 

rrii 
L 



Pxi i 



yi ^Pyi ' Pzi F z - 



P, 



P, = 



N pi 



gk B T 



(146) 
(147) 

(148) 

(149) 
(150) 
(151) 
(152) 

(153) 
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f Pl 



Pl = sA 



1 / N n 2 N \ 



K i=l ' i=i 

The Hamiltonian in Eq. (|145fl is separated into six terms as follows: 



(154) 
(155) 



#NPA 



1-RB1 



N ~9 N „/2 i 12 

2miS 2 L 2 ^ 2m l s 2 

i—l i—l 
N 



i=l ^ ' 



Hnpa 



1-RB2 



^NPA1-RB3 

-HnPAI — RB4 
-HnPA1-RB5 

#NPA1-RB6 



i=l 
JV 

8/ \ " 



p2 

'2Q ' 



(156) 
(157) 

(158) 

(159) 
(160) 

(161) 



In order to obtain the second-order symplectic formula, the time propagator exp [DAt] is again decomposed to a 
product of six time propagators as in Eq. (|106f) . The symplectic time developments are then given by 



exp [DxAt] 



TTLiSL 2 
,V2 



At, 



eatp [DiA*] i/j = Vl + ^-At 

mis 

J2 



exp [D 1 At] Zi = Zi + -^-At , 
rriiS 

exp[DiAt]qj = cos (CiiAt) g 4 + sin (Cii At) "Pi*^ , where £a 



4fis 



exp [DiAt] ^ = cos (CaAi) 7Tj + sin (C;iAi) Pitt- 

fJL f? N ,2 ^ 

[Cl A t ]P, = P. + S-^fe+VJ^ 



exp 



2m,iS 2 L 2 ^ 2rriiS 2 

\i—l i—l 
N 



■ ^ 2/iC?! - 5 feTo log s + H - gk B T At 



N 



^2 



exp [A At] P L = P L + V At , 



(162) 
(163) 
(164) 
(165) 
(166) 

(167) 
(168) 
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exp [D 2 At] q, = cos (C i2 A*) q t + sin {( l2 At) V 2 q l , where ( a = j^^Vtqi , (169) 

exp [D 2 At] Tt\ = cos(Ci 2 At)7r- +sin(Ci2At)P 2 7r- , (170) 

exp [D 2 At] P s = P s + ^ 27 2 C? 2 J At , (171) 

exp [D 3 At] q, = cos (( t3 At) q, + sin (Q 3 At) V 3 q t , where ( l3 = — Trf P 3 9i , (172) 

47 3 s 

exp [,D 3 At] < = cos (C, 3 At) tt< + sin (C l3 A£) V 3 n' t , (173) 

exp [D 3 At] P s = P s + 2I 3(k\ At > ( 174 ) 

exp [D 4 At] P s - P s - ||Ai , (175) 

exp [D A At] L = L + s^At , (176) 

exp [D 5 At] p xi = p xl + sLF xl At, (177) 

exp [D a At]p' yi = p' yi + sF yi At , (178) 

exp [D 5 At]p' zi = p' zi + sF zi At , (179) 

exp [D 5 At] iz't = Tv' l + 2sS(q l )N\ 4) At , (180) 

exp [D 5 At] P s = P S -{E + P Q AL) At , (181) 



\ i=l 



<>xp[D 5 At]P L = P h - ,s ^- > ^ /•;. - P U A j A/ . (1S2) 

ii 

2Q 

exp [D 6 At] F s = P s / ( 1 + -^At ) . (184) 



exp [D 6 A<] s = s [ 1 + -2-At ) , (183) 



2g 



F. Symplectic molecular dynamics algorithm for rigid-body molecules combined with the Nose-Poincare 

thermostat and the Parrinello-Rahman barostat 

In this subsection we present an explicit symplectic MD algorithm for rigid-body molecules in the isothermal-isobaric 
ensemble with simulation-cell deformation. The Hamiltonian is given by combing the Hamiltonian in Eq. I|23|) and 
the Parrinello-Rahman barostat [lj| as follows: 



H- 



NPPR-RB 



.1=1 
pi 



N 1 „-l N 1 „ „ 



(185) 



where L is the matrix of cell parameters, Pl is the conjugate momenta for L, and G is given by L L . The scaled 
momentum p i and the scaled coordinate fi are related to p i and r% here by 
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The equations of motion are given by 



Pi 



= Fi-- Pi -[ LL 



(4) 



P, 



Pi 



Q 



P s = E — + E - gknTo , 



L 

Pl 



S 

W 



N 



N 



(188) 

(189) 
(190) 
(191) 
(192) 

(193) 

(194) 

(195) 



where er is related to L by er = VL and 1 is the identity matrix. Note that p iP J and FirJ are dyadic tensors, 
whose {a, (3) elements (a,/3 = x,y,z) are p a iP/3i and F a irpi, respectively. The Hamiltonian in Eq. (|185|) is also 
separated into six terms as follows: 



H- 



NPPR-RBl 



' N 1 N 1 

E^^ G p*+E 



iV 



Irrns 1 
1 



8/is 2 



Trf I + gknTo log s - H 



NPPR-RB2 



8I 2 S 2 \ 1 



^NPPR— RB3 

-HnPPR— RB4 
-£?NPPR-RB5 
^NPPR— RB6 



8I 3 s 2 \ 1 



•E 

iV 

*E 

i=l 

— Tr I P T P L 

2W \ L 

s\E(r( N \q( N \L) + P Q V 

p2 

a-2- . 

2g 



(196) 
(197) 

(198) 

(199) 
(200) 
(201) 
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The symplectic time developments are given using the decomposition of exp [P/At] in Eq. i|106|) by 
exp [DiAi] fj 



. At-- 1 . 
r% H G Pi , 



exp [D 1 At] q % = cos (Cu At) q % + sin (Ci At) Vrfi , where Qi 
exp [P^At] tt^ = cos (Ca At) ttJ + sin (C«At) Pi< , 



4/i s 



JV 



1 



exp^At]^ = P S +(£__^ G - 



2V \ 

5^ 2/xC?! - gk^Ta log s + P - fffeTo At 
i=i / 



exp [£>! At] P L = 



r> At I ^ , 



1 

TO;S 



Pi 



Pi 



cos (C i2 At) + sin (Ci2 At) P 2 q l , where Ci2 = -— Tr- 1 ^^ , 

4i 2 s 



P s + ^ 2I 2 C? 2 j At 



cos (C i3 At) q % + sin {( i3 At) V 3 q, , where Ci3 = — < T P 3 g 2 , 

4/3 s 



exp [D 2 At] q, 

exp [P> 2 At] tt^ = cos(Ci 2 At)7r^+sm(Ci2At)'P27r- 
exp [D 2 At] P s 

exp [D 3 At] q t 
exp [P> 3 At] Tv'i = cos At) tv[ + sin (foAt) P 3 < , 
exp [D 3 At] P s = P s + 2 hC^j A* , 

At /~ T - 
exp [£» 4 At] P s = P s - — Tr f P L P L 

<-> <-> s At — 
exp [D 4 At] i = L + — P L , 

cxp[D 5 At]p l = p z + sL FiAt, 
exp [£> 5 At] Tr'i = 7v' l + 2sS(q l )N ( i i) At 



exp[D a At]P 3 = P s 



E(rl N \qW,L)+P V 



At 



N 



exp [D 5 At] P L = P L +sAt ( - FirJ - P 1 



P 



; = 1 
2 



exp [D 6 At] s = s ( 1 + gA/ 
<'-M>[A : A/]P. : P,/ ( l + ||At 



(202) 
(203) 
(204) 

(205) 

(206) 

(207) 
(208) 
(209) 

(210) 
(211) 
(212) 

(213) 

(214) 

(215) 
(216) 
(217) 

(218) 

(219) 
(220) 



G. Symplectic condition and time reversibility 

In this section we discuss the symplectic condition and the time reversibility |9J . Let us suppose a time- independent 
canonical transformation from 



(221) 
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to 

r--(2g:3). « 

where Q and P arc the generalized coordinate and the generalized momentum, respectively. The canonical equation 
of r is given by 

dH , . 

r = J w , (223) 



where 



1 
1 



(224) 



Because r is given by the canonical transformation from F, the canonical equation of T' is also given by 



r' = J§p • (225) 



The time derivative of r'(T) is derived in another way by the chain rule: 

r' = |r = Mr = mj^ = mjm^ , (226) 

where M is the Jacobian matrix for the canonical transformation from T to T' and its («, j) element is given by 

dr' 



Mij = ^ . (227) 



or, 



Comparing Eq. (|225ll and Eq. I|22t)|) . we obtain the symplectic condition: 



MJM T = J . (228) 

In general, the generalized coordinates and momenta obtained by a Hamiltonian dynamics fulfills the symplectic 
condition in Eq. (|228|l . 

Each factor in the decompositions in Eqs. I|45[) and (|106fl is a time propagator based on the corresponding Hamilto- 
nian. For example, exp [DiAi] in Eqs. Ij45(l is a time propagator by the Hamiltonian of Hnp-rbi. Therefore, the time 
developments by the decompositions in Eqs. (|45|l and l|lUtj|) fulfill the symplectic condition. All variables in Eqs. i|56|) - 
(I8GI1 are canonical variables such as r;, p' i} q i , ir^, s, and P s . Besides, the time propagator here is decomposed so that 
the MD algorithm will be time reversible, namely, exp [— DAt] exp [£>Ai] = 1 holds in Eqs. H45|) and (|106fl . 

Employing the symplectic MD algorithm, there is a conserved quantity which is close to the Hamiltonian Q. It 
means that the long-time deviation of the Hamiltonian is suppressed. Therefore, we can perform a MD simulation 
more stably than by conventional nonsymplectic algorithms. 

From the symplectic condition in Eq. I|228|) . the Jacobian determinant is calculated as one: 

detM = 1 . (229) 

It means that the phase-space volume is conserved during the simulation. Note that the phase-space-volume con- 
servation is a necessary condition of the symplectic condition and not a sufficient condition. The condition that the 
Jacobian determinant is one does not always mean symplectic. Even if the Jacobian determinant is one, there is 
not always a conserved quantity which is close to the Hamiltonian. In other words, there are nonsymplectic MD 
algorithms which are phase-space volume conserving and time reversible. The time propagators in these nonsymplec- 
tic algorithms are not based on Hamiltonian and the variables are not canonical variables. That is, the symplectic 
condition in Eq. (|228|) is not fulfilled. Therefore, there is not a conserved quantity which is close to the Hamiltonian. 
It means that the value of the Hamiltonian deviates gradually from its initial value in a long-time simulation. In the 
next section we compare our symplectic algorithm with the nonsymplectic time-reversible algorithms. 



III. COMPARISONS WITH NONSYMPLECTIC TIME-REVERSIBLE ALGORITHMS 

In this section we explain three nonsymplectic algorithms in the canonical ensemble, which are time reversible. We 
then apply our symplectic algorithm and these nonsymplectic algorithms to a rigid-body water model and compare 
them numerically. 
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A. Molecular dynamics algorithm based on the Nose-Poincare Thermostat and the Nonsymplectic 

Rigid-Body Algorithm 



Instead of the symplectic rigid-body MD algorithm by Miller et al. @, we here combine the nonsymplectic rigid- 
body MD algorithm by Matubayasi and Nakahara Q with the Nose-Poincare thermostat [TlUT^. In this algorithm, 
angular velocity u)^ = suji instead of 7r£ is employed, that is, the variables here are n, p' i: q i , w'-, s, and P s . 

The time propagator exp [DAt] is decomposed by 
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where each time propagator is given by 
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(230) 



(231) 

(232) 

(233) 

(234) 
(235) 



B. Molecular dynamics algorithm based on the Nose-Hoover Thermostat and the Symplectic Rigid-Body 

Algorithm 



We here combine the symplectic rigid-body MD algorithm [f| with the Nose-Hoover thermostat 0, S EH (the 
latter is nonsymplectic). This combination has been employed in Ref. |15|| . Instead of s and P s , r) = logs and 
£ = P s /Q are used for the thermostat, that is, the variables employed here are r^, p i: q i7 rj, and £. 

The time propagator exp [DAt] is decomposed by [ljj 
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(236) 
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where each time propagator is given by [if 
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, mi dri 
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9^ 



(237) 



* - £ 35 (55*) ■ 4 + g i (^,) ■ a , (238) 

W N 

N „ AT „ „ 



^i ^ dr] 



1 / ^ ? N 1 T \ (9 

D « = q (£ ^ + £ inJSiqjDiS (g,)w< - .gfeToj - . (242) 

We remark that we can also make another second-order integrator by the decomposition in Eq. p()6|l instead of 
Eq. l|286[l . However, the original time reversible algorithm for the Nose-Hoover thermostat decomposed the time 
propagator as in Eq. i|23(i[) |l0|. thus we used this decomposition. 



C. Molecular dynamics algorithm based on the Nose-Hoover Thermostat and the Nonsymplectic 

Rigid-Body Algorithm 

We can also make a nonsymplectic algorithm by the rigid-body algorithm by Matubayasi and Nakahara and the 
Nose- Hoover thermostat 0, & • I n this algorithm the following variables are developed with time: p i: q t , 
Ui, r), and £. 

The time propagator exp [DAt] is decomposed as in Eq. I|236|) . Each decomposed time propagator is given by 

ft = E£-£+EU^>rV£. p«) 
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D. Numerical comparisons: application to a pure water system 

We applied the symplectic and nonsymplectic MD algorithms to a rigid-body model of water in the canonical 
ensemble. We employed the TIP3P rigid-body model for the water molecules We used 80 water molecules in 
a cubic unit cell with periodic boundary conditions. The temperature was set at 300 K and the mass density was 
set to 0.997 g/cm 3 . The electrostatic potential was calculated by the Ewald method. We calculated the van der 
Waals interaction, which is given by the Lennerd- Jones term, of all pairs of the molecules within the minimum image 
convention instead of introducing the spherical potential cutoff. We tested the time steps of At — 2 fs, 3 fs, 4 fs, 
and 5 fs. We performed the MD simulations for 1.5 ns in all cases of At. We employed Eas. l|56|) - (|86[) for the Nose- 
Poincare thermostat and symplectic rigid-body MD simulations, Eas. l231|l - (|235fl for the Nose-Poincare thermostat 
and nonsymplectic rigid-body MD simulations, Eas. (|237(l - Ij242(l for the Nose-Hoover thermostat and symplectic 
rigid-body MD simulations, and the time development in Eas. l|243(l - Ij248(l for the Nose-Hoover thermostat and 
nonsymplectic rigid-body MD simulations. The same initial conditions were used for all algorithms and time steps. 

We observed the deviations of the Nose Hamiltonian from its initial values: 

N „/2 N 1 «_ „ p2 

6H ® = E 2^ + £ g^^^DiS (gM + £(rW,?W) + ^ + gk B T logs - H . (249) 

Figures H 12 El and^Jshow SH(t) for At — 2 fs, 3 fs, 4 fs, and 5 fs, respectively. The gradient of the linear fitting for 
each 8H(t) is shown in Table [I] 

In every nonsymplectic MD algorithm, Hamiltonian deviates from its initial value as time passes even for At = 2 f s 
as shown in Figs. IHb)-(d). This deviation increases as the time step increases from At = 2 fs to 5 fs. Note that the 
energy scale in the ordinate increases as At increases. 

On the other hand, the Nose-Poincare thermostat and symplectic rigid-body MD algorithm guarantees the existence 
of a conserved quantity which is close to the Hamiltonian. Because of this conserved quantity, the Hamiltonian was 
conserved well for time steps of At = 2 fs, 3 fs, and 4 fs as shown in Figs. [IJ a )"ISt a )' The Hamiltonian starts to 
deviate slightly by dSH(t)/dt = 3.7 x 10~ 3 kcal/mol/ns in the case of At = 5 fs as shown in Tabled However, it is 
only two or five percent of SH for the other nonsymplectic MD simulations (TableOJ. This fact means that employing 
the combination of the Nose-Poincare thermostat and the symplectic rigid-body algorithm, one can take a time step 
of as much as 4 fs. This time step is longer than typical values of 0.5 fs to 2 fs which are used by the conventional 
nonsymplectic algorithms. 

IV. CONCLUSIONS 

We have proposed an explicit symplectic MD algorithm for rigid-body molecules in the canonical ensemble. This 
algorithm is based on the Nose-Poincare thermostat 0, ^| and the symplectic rigid-body algorithm [f| . We also 
have presented an explicit symplectic MD algorithm for rigid-body molecules in the isothermal-isobaric ensembles by 
combining the Andersen barostat |16| with the symplectic algorithm in the canonical ensemble. As a modification of 
the isothermal-isobaric algorithm, we further presented the symplectic integrator in the constant normal pressure and 
lateral surface area ensemble and a symplectic algorithm combined with the Parrinello- Rahman algorithm. Employing 
the symplectic MD algorithm, there is a conserved quantity which is close to the Hamiltonian. Therefore, we can 
perform a MD simulation more stably than by conventional nonsymplectic algorithms. 

In order to establish this fact numerically, we have applied this algorithm to a TIP3P pure water system at 300 K and 
compared the time evolution of the Hamiltonian with those by the nonsymplectic algorithms. These nonsymplectic 
algorithms are based on the Nose-Poincare thermostat [ill Il2| | and the nonsymplectic rigid-body algorithm Q , based 
on the Nose-Hoover thermostat ^(j an d the symplectic rigid-body algorithm and based on the Nose-Hoover 
thermostat ^(j an d the nonsymplectic rigid-body algorithm Q. In these nonsymplectic algorithms, the Hamiltonian 
deviates gradually from its initial value in all cases of the time steps At = 2 fs, 3 fs, 4 fs, and 5 fs. On the other hand, 
the Hamiltonian was conserved well even for a time step of 4 fs in our symplectic algorithm. 

The rigid-body model for molecules can be employed not only for a water system but also for a biomolecular system. 
For example, a partial rigid-body model |l5j is often used for a part of a peptide and a protein, in particular, for a 
hydrogen-including part to alleviate a fast motion of the hydrogen atom. Our algorithms will thus be of great use for 
MD simulations of an aqueous solution and a biomolecular system at a constant temperature and/or pressure. 
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TABLE I: Drift of the Hamiltonian per nanosecond dSH/dt (kcal/mol/ns). 
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3 fs 




4 fs 




5 fs 


Nose-Poincare and symplectic rigid-body MD 


-1.4 x 10" 
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Nose-Poincare and nonsymplectic rigid-body MD 
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Nose-Hoover and symplectic rigid-body MD 
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Nose-Hoover and nonsymplectic rigid-body MD 
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FIG. 1: The time series of the difference of Hamiltonian from its initial value 8H(t). The time step was set to At = 2 fs. (a) 
Nose-Poincare thermostat and symplectic rigid-body MD, (b) Nose-Poincare thermostat and nonsymplectic rigid-body MD, (c) 
Nose-Hoover thermostat and symplectic rigid-body MD, and (d) Nose-Hoover thermostat and nonsymplectic rigid-body MD. 
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FIG. 2: The time series of 5H(t). The time step was set to At — 3 fs. See the caption of Fig. 1 for further details. 
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FIG. 3: The time series of 8H(t). The time step was set to At — 4 fs. See the caption of Fig. 1 for further details 
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FIG. 4: The time series of 5H(t). The time step was set to At — 5 fs. See the caption of Fig. 1 for further details 



